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The effect of electronic current on the atomic motion still poses many open questions, and several 
mechanisms are at play. Recently there has been focus on the importance of the current-induced non- 
conservative forces (NC) and Berry-phase derived forces (BP) regarding the stability of molecular- 
scale contacts. Systems based on molecules bridging electrically gated graphene electrodes may 
offer an interesting test-bed for these effects. We employ a semi-classical Langevin approach in 
combination with DFT calculations to study the current-induced vibrational dynamics of an atomic 
carbon chain connecting electrically gated graphene electrodes. This illustrates how the device 
stability can be predicted solely from the modes obtained from the Langevin equation including 
the current induced forces. We point out that the gate offers control of the current independent of 
bias voltage which can be used to explore current-induced vibrational instabilities due the NC/BP 
forces. Furthermore, using tight-binding and the Brenner potential we illustrate how Langevin-type 
molecular dynamics can be performed including the Joule heating effect for the carbon chain systems. 
Molecular dynamics including current-induced forces enables an energy redistribution mechanism 
among the modes, mediated by anharmonic interactions, which is found to be vital in the description 
of the electronic heating. We have developed a semi-classical Langevin equation approach which can 
be used to explore current-induced dynamics and instabilities. We find instabilities at experimentally 
relevant bias and gate voltages for the carbon chain system. 

PACS numbers: 



I. INTRODUCTION 

The consequences of electronic current on the motion of atoms have become accentuated with the on- 
going quest for molecular-scale electronics [DHJ- The atomic motion due to electrical current is behind 
the long-term breakdown of interconnects leading to failure in integrated circuits. This effect is of even 
greater importance for systems where the bottle-neck for the electronic current is a few chemical bonds. 
The inelastic scattering by electrons on atomic vibrations leads to the well-known Joule heating, which 
can have impact on the electronic behavior and stability. However, recently it has been pointed out [5 
[H] that other current-induced forces can play a role. For instance, in the case of molecular contacts 
with conductance on the order of Go = 2e 2 /h = l/12.9kSl (e being the electron charge and h Planck's 
constant), and under "high" bias voltage (~ IV), the current-induced forces which does not conserve the 
energy of the atomic motion may lead to run-away behavior. However, experiments in this regime are 
very challenging. For example, for the typical experiments involving molecular-scale contacts between 
bulk electrodes it is not possible to image the atomic structure in the presence of contact and current. 
Furthermore it is highly non-trivial to add additional gate potentials in order to modify the electronic 
structure and get independent control of bias voltage and current [21 [H] ■ 

On the theoretical side, it is desirable to develop computer simulations techniques such as molecular 
dynamics (MD), preferably without adjustable parameters, to study in detail the current-driven complex 
atomic processes. To this end, we have recently developed an approach based on the semi-classical 
Langevin equation, which may form the basis of MD. In this description the non-equilibrium electronic 
environment is described like an effective "bath" influencing the atomic dynamics. Especially, we identify 
the forces acting on the atoms due to the electronic current. These include "extra" fluctuating forces 
yielding the Joule heating, a non-conservative " electron- wind" force (denoted NC), recently discussed 
by Todorov and co-workers [5], and a Lorentz-like force originating from the quantum- mechanical "Berry 
phase" of the electronic subsystem^ (denoted BP). The purpose of this article is two-fold. We will 
illustrate this semi-classical Langevin approach, and show how the current-induced effects could be 
investigated in molecular contacts connecting gated graphene or nanotube electrodes. 

Graphene is now being explored very extensively due to its eminent electrical and thermal transport 
properties [TUHH]. Besides being highly important in their own right, carbon nanotube or graphene- 
based nanostructures may offer an interesting test bed for studies of current-induced effects at the 
atomic scale. For such systems experiments with atomic resolution, using for instance state-of-the-art 
electron microscopes, can be performed in the presence of current, allowing the dynamics to be followed 



down to single adatoms[T3]- Electronic current has been used to induce changes in graphene-edges, which 
were monitored while simultaneously passing a current through the structure [14j. This was explained as 
carbon edge-dimers desorbing due to Joule-heating [T5] . Taking this a step further one can imagine that 
nano-structured nanotubes or graphene can be used as an electrode interface to molecular devices based 
on organic chemistry [T5]. Especially promising prospects include the inherent 2D geometry of graphene 
which both enables straight-forward electrical gating, and atomic-scale imaging in the presence of current. 
There has been a number of microscopy studies of single-atom carbon chains bridging graphene |13l [T7] 
or nanotubes [TB]- On the theoretical side various aspects of these systems have been studied such as 
the formation of chains [HI [20], their stability [2 1|. and electron transport properties [22"k24| . Here we 
will explore the current-induced forces and nano-scale Joule heating of the carbon chain system between 
electrically gated graphene electrodes. 

The paper is organized as follows. After a brief sketch of the semi-classical Langevin method, we will 
use it to study the dynamics of the carbon chain as a function of bias and gate voltages. We point out 
that the gate, which offers independent control of bias voltage and current in the system, can be used 
to explore current-induced vibrational instabilities in the current-carrying chain. Finally, we illustrate 
how the Langevin molecular dynamics can be performed including the Joule heating effect for a carbon 
chain system, using tight-binding and the Brenner potential. 

II. SEMI-CLASSICAL LANGEVIN DYNAMICS 

We will here sketch the Langevin approach. For a classical oscillator system (mass-scaled coordinate 
x) in a general non-linear force-field, F, coupled linearly to a bath of harmonic oscillators it is possible 
to eliminate the bath variables and describe it using the generalized Langevin equation, [2"51 - f2"T] . 

x = F(x)- f W{t,t')x(t')dt' + £{t). (1) 

Here the bath is influencing the motion via two distinct force contributions, (i) a retarded time-kernel, 
n r , describing the back-action at time t after propagation in the bath due to the motion of x at previous 
times, and (ii) a force £ of statistical nature originating from the thermal fluctuations of the bath. In 
the case of thermodynamic equilibrium £ is characterized by a temperature and is related to n r by the 
fluctuation-dissipation theorem. Note that in general x, F, and £ are vectors and n r is a matrix. This 
method has been used by Wang and co-workers |28, [29] to describe thermal transport in the quantum 
limit, using phonons in the two connecting reservoirs with different temperature as baths and by including 
their quantum fluctuations in £. This reproduce the Landauer result of thermal transport in the harmonic 
case [33]. 

It is possible to reach a semi-classical Langevin equation description of the motion of the ions coupled 
to the electron gas if we assume a coupling to the electronic environment which is linear - either in the 
displacement from an equilibrium or velocity (adiabatic expansion) of the ions. This Langevin/Brownian 
motion approach to atom scattering at metal surfaces has a rather long history in the case of metal 
electrons in thermal equilibrium [30l l3"T] . 

We have extended this to describe the dynamics of the ions in a nano-conductor between metal 
electrodes in the non-equilibrium case, where an electronic current is present [51 132). In order to sketch 
the derivation, we consider a displacement dependent tight-binding model with electron states in the 
scattering region of interest k, I, and H e i being the static electronic Hamiltonian (scattering region and 
its coupling to left and right electrodes [3"5] ). 

H = H ph {x) +H el + ^2 M nM clci x n . (2) 

Here a; is a column vector made of mass-normalised displacement operator of each degrees of freedom, e.g., 
x n = y/m n u n , u n and m n being the displacement operator and mass, H p h — |i x + \x T Kx corresponds 
to a set of harmonic oscillators that couple with the electrons, K being the dynamical matrix. We imagine 
a localized basis-set describing the electrons in the scattering region, where ct (c&) is the electron creation 
(annihilation) operator at site k in this region [H3]. Here we only consider the coupling to the electronic 
bath, but the linear coupling to an external phonon bath can be taken into account along the same lines 
and adds a contribution to n r . The derivation and result for linearly coupled harmonic phonon bath 
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is similar, and has been given in Ref. [25]. Alternatively, the dynamics of some external phonons, not 
coupling to the electrons directly, may be treated explicitly in actual MD calculations, as we illustrate 
later (regions DL,DR in Fig. [6^,). The electron-phonon coupling corresponds to matrix elements of the 
force operator M ni jy = (k\S/ Xn H e i\l) . We have assumed that M is small by keeping only the term linear 
in x. 

We may obtain an equation of motion for x using Heisenberg's equation of motion, x = i[x, H], using 
atomic units (h = 1) and implicit mode index (n), 

x = -Kx - Mfci4ci = -Kx + f e . (3) 

kl 

The latter term, f e , describes the "forces" due to the interaction with the electron gas. Importantly, 
these forces have a random nature[35]. We can calculate the mean value of f e by averaging it over the 
non-equilibrium electronic state, 

(f e ) = zTr[MG<(M)] = Tr[(-V x H e i) p(t)} . (4) 

Here we have introduced the electronic lesser-Greens function, G^(i, t) = i(c](t)cj(t)) , which is equiva- 
lent to the density matrix, p (times —i), and depends on x(t), since the electrons are coupled to the x in 
the Hamiltonian. This is similar to the expression for an average force in Ehrenfest dynamics [5]. 

We can evaluate this perturbatively using the unperturbed electron lesser Green's function, G^ , cor- 
responding to the case of steady-state electron transport without electron-phonon interaction[33 , 

Go M = iA L (uj)n F (uj - p L ) + iA R (uj)n F (uj - p R ) , (5) 

where A F / R are density of state matrices for electronic states originating in the left/right electrodes, 
each with chemical potential pr /t?[33]. which differ for finite bias voltage, V, as ph — PR = eV , and 
n F (uj) — l/(e u ' kBT + 1) is the Fermi-Dirac distribution function. We thus treat the non-equilibrium 
electronic system as a reservoir unperturbed by the phonons. Using the non-equilibrium Greens function 
(NEGF) technique [5^]. we may write the 2nd lowest orders in M of (f e ) as, 

(f e (t)) w J duTr[(-V x H el )p ] - jf U r (t,t')x(t')dt' . (6) 

The first term yields a constant force due to the change in electronic bonding with bias and a "direct 
force" due to interaction of charges with the field [37] . Here po — p cq + 5p is the non-equilibrium electron 
density matrix without electron-phonon interaction. We split it into an equilibrium contribution p eq and 
a nonequilibrium correction 8 p. In linear response, we get a term E ■ x from the field in H e i, £ being 
the external field, yielding a "direct" force involving the equilibrium p cq . We also get a term involving 
H e i {£ — 0) together with the change in density to first order in the field Ap oc £ in the first term of ([6| 
resulting from the change of density in the chemical bonds due to the current [38, 39J. 

The second contribution is the retarded back-action of the electron gas due to the motion and is 
equivalent to the retarded phonon self-energy. In the steady state, n r only depends on the time difference, 
and it is convenient to work in the frequency(energy) domain. This can be expressed using the coupling- 
weighted electron-hole pair density of states, A"' 9 , inside or between electrodes a, (3 £ L, R, 

n r (t-t') = J n r ( u ) e -^, (7) 
n»= / , A(a/) .M , (8) 

J ur — uj — id 

where A can be expressed in terms of the electrode DOS, 

A = 53 A°^, (9) 

a/3 

A°i(w) = / — Tr[M m A a (o;')M„A^a;' - uj)} {n F {J - p a ) - n F {J - uj - p fj )) . (10) 



We have included a factor of 2 from the spin-degeneracy and have explicitly included the mode index, 
to, n on the coupling matrices, M, and on A in Eq. |10[ 
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The forces described by n^„(a;) in Eq. [6] contains a number of interesting current- induced effects, ft 
is instructive to split the kernel into parts, 

U r mn (uj) = wrRe(A mn (u;)) - 7rIm(A m „(w)) + nH{Re(A mn )}(uj) + i7rK{Im(A mn )}(o;), (f f ) 

where T-L{f{x')}{x) — J g J, x _l dx' is the Hilbert transform. The A matrix has the following symmetry 
properties when exchanging modes(n O m) and electrodes( a O /3), 

a^h = a^;h, (12) 

and 

A£f» - -Af« (-«), (13) 



which arc helpful when examining the terms in Eq. ( 11 ), which are summarized in the following: Friction 



The first term in Eq. 11 is imaginary and symmetric in mode index m,n. It describes the friction force, 
due to the generation of electron-hole pairs in the electronic environment by the ionic motion. This 
process exists even in equilibrium 31 . For slowly varying A. L / R with energy compared to the vibrational 
energies (wide-band limit) we obtain the simple time-local electronic friction force, —rj e ix, with 

^ = _ttRc(A) 1^ 

U) Z7T 

a,0 



NC (wind) force- The second term in Eq. 11 is real and anti-symmetric, which means that the general 
curl of this force is not zero. It is describing the NC force, discussed very recently by Dundas and co- 
workers [5^. This force is finite, even in the limit of zero frequency, where the friction and Joule heating 
effect is not important anymore. 

Renormalization-The third term is real and symmetric and can be interpreted as a renormalization 
of the dynamical matrix. It contains an equilibrium part and a nonequilibrium correction. The equilib- 
rium part is already included in the dynamical matrix if wc calculate it within the Born-Oppcnhcimer 
approximation. The nonequilibrium part gives a bias-induced modification of the harmonic potential. 

BP force -Finally, the last term is imaginary, anti-symmetric, and proportional to tu for small frequen- 
cies. It can be identified as the "Berry phase" (BP) force in Ref. [5]. Since the direction of this force is 
always normal to the velocity in the abstract phase space, it does no work resembling a Lorentz force 
with effective magnetic field 

B=/ WA W) }M (15) 

Random forces- The randomness of the force f e is characterized by its correlation function in fre- 
quency domain, which can again be calculated using NEGF. However, we note that since f e is a quan- 
tum operator, (/ e (i)/ e (0)) does not result in a real number. Instead we use the symmetrized and real 
(fe(t)fe(0) + (/ e (0)/J(i)) T ). This expression equals the semi-classical result obtained from the path- 
integral derivation of the Langevin equationJB, 35] and reads in Fourier space, 

= </e/ e T >M - </ e )</e) T H = -TT^COth ( CJ "^;~^ ) ) A«» . (16) 

This spectral power density can be used to generate an instance of the Gaussian random noise as a 
function of time which is needed in MD simulations. Most importantly this random force contains 
not only the thermal excitations but also the excess excitations leading to Joule heating [5^] via the 
dependence of the chemical potentials ^l—^r = eV . Thus with this formalism it is possible to disentangle 
the various contributions to the forces, being either deterministic or random in nature. 

III. CURRENT-INDUCED VIBRATIONAL INSTABILITY 

We now turn to illustrations of the use of the semi-classical Langevin equation to describe current- 
induced effects. In this section we employ it to study the effect of the current-induced forces and 
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FIG. 1: The system considered in the present study is a four-atom carbon chain bridging two graphene electrodes. 
The dangling bonds are passivated by Hydrogen atoms. In addition to the bias applied between the left (L) and 
right (R) electrodes (Vb), a gate potential (V^) can also be applied perpendicular to the graphene surface. The 
center panel shows the calculated contour plot of the electrostatic potential drop across the junction at V g = OV, 
and Vb = IV. The equal drop at the left and right electrodes reflects the electron-hole symmetry for V g = OV 40 . 
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FIG. 2: Current- Voltage (7 — Vb) curves at different V a . 



Joule heating on the stability of the system within the harmonic approximation. We will here ignore the 
coupling to electrode phonons. This makes an eigen-mode analysis possible which ease the interpretation 
of the results. The model system we use is shown in Fig. [T] where a four-atom carbon chain is bridged 
between two graphene electrodes(L and R). We assume a field effect transistor setup, where a gate 
potential, V g , is applied to the system in addition to the bias applied between the two electrodes, Vb- We 
will show that this offers a convenient way to explore current-induced vibrational instabilities. We can 
already see the effect of the gate potential in the current-voltage (I — Vb) characteristics shown in Fig. [2] 
The effect of the NC and BP forces is to couple different phonon modes with nearly similar frequen- 
cies. From now on, we will focus on the two phonon modes around 200 meV, shown in Fig. [3j since 
the alternating-bond-length-type modes (200 meV) couples most strongly with the electrical current. 
This type of modes also gives rise to the most intensive Raman signals in unpassivated chains between 
graphene- like pieces [H]. 

The calculation is performed using the SIESTA density-functional theory (DFT) method [42], which 
has been extended to study elastic j55] and inelastic [M] transport in molecular conductors. We use 
similar parameters as detailed in Ref. [33], and in order to keep the calculation simple and tractable, we 
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FIG. 3: (a) Motion of the two phonon modes around 200 meV. (b) Motion of the runaway mode at V g = 0.6V, 
and Vb = IV. We depict the motion using a number of discrete time steps roughly corresponding to a full period. 
The position of each atom is depicted as a circle for a sequence of time steps indicated by an increasing radius with 
time. The motion is a phase-shifted linear combination of the two modes in (a). We can see the elliptical motion 
of the carbon atoms from the plot. The inclosed area indicates that work can be done by the current-induced 
NC force. 





FIG. 4: (a) Inverse Q- factor (1/Q) as a function of gate voltage, V g , at Vb = IV for the two modes around 200 
meV. (b) 1/Q as a function of bias voltage, Vb, at fixed gate voltage V g = 0.6V, for the same pair of phonon 
modes. 








FIG. 5: (a) Effective phonon number (N) for the two phonon modes around 200meV as a function of gate voltage, 
V g , at fixed bias voltage, Vb = IV. (b) N as a function of bias voltage, 14, at fixed gate voltage V g — 0.6V. Note 
that it diverges at the critical point when the damping (1/Q) in Fig. [4] goes to zero. 



model the electrodes by just employing the T k-point in the transverse electrode direction. The electron- 
phonon coupling matrix (M) is calculated at zero bias, while we calculate the electronic structure at 
finite bias. We note that the voltage dependence of the coupling matrix could play a role, but this is 
beyond the present more illustrative purpose [43j. Based on these approximations, we can calculate the 
full w-dependent A function, and the self-energies, II r . To perform the eigen-mode analysis, we further 
assume linear w-dependent friction, Berry force (BP), constant non-conservative force (NC), and ignore 
the renormalization of the dynamical matrix. 

We model the effect of 14 as a shift of the equilibrium chemical potential, Ep. In this way we can 
tune the electronic structure within the bias window by changing the gate potential. In the following, 
we will look at the bias and gate dependence of the inverse Q-factor (1/Q) and effective phonon number 
N . The inverse Q-factor for mode i (note we use index i for full modes including the current-induced 
forces) is defined as 

(17) 

Re-fa;,} 

where u)i is the eigenvalues of the full dynamical matrix, including the current-induced forces. These 
modes will then consist of linear combinations of the "unperturbed" normal modes of the system, n, m, as 
calculated using the standard Born-Oppenheimer approximation. The phonon number can be calculated 
from the displacement correlation function, 

1 f doj 

N i + 2~ Re ^ / ( X ^>H^:- ( 18 ) 

We show the bias and gate potential dependence of the inverse Q-factor and phonon number in Fig. [4] 
and Fig. [5] The coupling of these two modes due to the bias (gate) dependent NC and BP force changes 
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their lifetime. The two modes always have opposite dependence. The vibrational instability happens at 
the critical point where 1/Q = around V g — ±0.4V. This corresponds to an infinite phonon number in 
Fig. [5j and we therefore call it a "runaway" mode. The motion of this mode at VJ, = IV, V g = 0.6V is 
plotted in Fig.^b). We can observe the elliptical motion in real-space of several atoms. This is critical 
because in order for the non-conservative force to do work on the atoms their motion has to enclose a 
finite area, either in real or in abstract phase space. 

Finally, we should mention that once hitting the instability threshold the current will drive the system 
to some highly anharmonic regime, where present eigenanalysis breaks down. One scenario is that the 
motion of the system will reach a limit cycle determined by the detailed anharmonic potential and the 
interaction with the current [7]. In this regime the detailed damping due to the coupling with phonons in 
the electrodes could be important, as well as the electron-phonon coupling could change from the value 
around the harmonic equilibrium position. In order to address this regime we can perform molecular 
dynamics simulations, taking into account both the coupling between different modes and their coupling 
with the electrode phonons, in order to study how the system actually react due to the instability. 

IV. MOLECULAR DYNAMICS WITH JOULE HEATING 

Next we illustrate the use of the Langevin equation to perform molecular dynamics simulations of a 
carbon chain system in the presence of current in the simplest possible setting, but now including the 
coupling to electrode phonons. Therefore we abandon the DFT approach, and instead employ the widely 
used 7r-tight-binding model with hopping parameter f3 = 2.7eV, and the Brenner potential for calculations 
of the inter-atomic forces |33j. We consider the unpassivated structure in Fig. [6] The electron-phonon 
coupling is modeled by the Harrison scaling law|45j. (3 = 2.7eV(ao/d) 2 , determining how f3 is modified 
if the nearest neighbor distance, d, is changed from the equilibrium value, ao — 1.4A. The same model 
has recently been applied to study the effect of strain on the electronic structure of graphene[lS]. In the 
simulation model the coupling to electrode phonons by a friction parameter, r) p h, and a corresponding 
white equilibrium phonon noise (£, P h£,ph) = ^VphkBT on the L,R electrode regions. This is similar to 
the stochastic boundary conditions [27], where L, i?-atoms act as a boundary. The set up for the MD is 
shown in (Fig. [6^,). We include electrode regions without interaction with the current (DL, DR), and a 
device region (D) where the current density is highest and where the nonconservative forces and Joule 
heating is included. 

Furthermore, instead of using the full non-local time-kernel for the electrons in Eq. |14[ we use the 
wide-band approximation, and neglect the off-diagonal elements of the electron noise spectral power 
density, (£ e £j)(u;). The diagonal of the electron spectral power can be approximated by white noise 
in the high bias and wide-band limits, where variations in the electronic DOS is neglected [4 7 r The 
assumption of a white noise spectrum implies neglect of the equilibrium zero-point motion of the atoms, 
but most importantly here, it includes the Joule heating effects, 

(K.)n(WD(w) = 2(r )el ) nn k B T + Re(Tr[A L {fi L )M n A R ( f i R )M n ]) S (19) 

A factor of 2 should be included in the case of spin-degeneracy. Based on the velocity Verlet algorithm [48 j 
we have carried out MD simulations at a varying voltage bias for zero gate bias (V g = 0V), and phonon 
friction, r\ v h- The MD results are summarized in figure [Tjlb-f). We note that for the present system 
set-up the non-conservative force is found not to play a dominant role compared to the effect of Joule 
heating. The main insight we gain from the MD example here is that the anharmonic couplings are 
important and effective in redistributing the energy supplied by the non-equilibrium electrons. 

The approximate local phonon friction, 77^, can in general be expressed from the slope of the cor- 
responding phonon self-energy at zero frequency as for electrons, see Eq. |14| However, here we have 
simply varied its value around this in order to quantify the dependence of the local electronic heating 
in the device region on this parameter (Fig. [Fjja). The electronic heating of the chain is found not to 
depend much on the phonon friction when this is chosen sufficiently high. This is an appealing result, 
since it indicates that the electronic heating does not depend critically on the measurement setup, but 
mainly on the nature of the actual constriction. This seems to be true as long as the heat flow away 
from the contacts is sufficient to maintain the temperature of the heat baths, and the chain acts as a 
bottleneck for the heat conduction. However, we note that for heat conduction in the quantum limit it 
is important to go beyond the white band approximation and include realistic self-energies for the L, R 
electrode phonons |49j. This will be explored in future work. 
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FIG. 6: (a) Definition of system regions with different types of noise contributions. Leads (L, R) have a 
well-defined temperature specified from the phonon noise, the device (D) temperature will be defined from the 
electronic heating, and the intermediate regions (DL, DR) are free and will obtain a temperature from propagating 
noise. In the MD setup no atoms are held fixed but periodic boundary conditions are applied. The figure describes 
how the setup in which the local temperatures plotted in figure (c+e) should be understood, (b) Temperature of 
the regions as a function of phonon friction. (c,d) Obtained temperatures at different atoms within the harmonic 
approximation, (c) The simulations are run at T=300K and at eVJ, = leV, and (d) varying bias voltages. (e,f) 
Corresponding atomistic temperature distributions including the anharmonic interactions. The lead temperature 
can exceed the equilibrium bath temperature due to propagating noise. Especially the anharmonic interactions 
redistribute part of the energy from the modes in the chain to the bulk modes in the lead. 



Inspired by the equipartition theorem, we define a local temperature variable for the atoms (indexed 
by a) with mass, m ai 

T B (t) = -^(t)>. (20) 
3k B 

Comparing the obtained temperature distributions with(Fig. [6j:,d) and without(Fig. [6j2,f ) the anhar- 
monic interactions show that anharmonic couplings between the vibrational modes has a significant 
influence on the heat transport properties and local Joule heating of the system. The heating is less 
localized in the chain due to anharmonicity. This originates from the coupling of different modes and 
an increased coupling to the surroundings for configurations where the atoms are displaced from their 
equilibrium positions. Modes localized in the chain can be heated up to very high temperatures in the 
harmonic approximation. When anharmonic interactions are included the energy is redistributed and 
the modes are collectively heated up. 

The electron-phonon interaction is typically included through a Taylor expansion of the electronic 
Hamiltonian around the equilibrium positions (Eq. [2]) . Within the time-local white noise approximation 
it is possible to address the effect of change of electronic Hamiltonian and, especially electron-phonon 
coupling, on the motion. This amounts to updating the friction and noise on the fly along the path. This 
is possible for the simple parametrization used here. Our preliminary results based on this approximation 
show that the extra noise contribution from the higher-order couplings may significantly influence the 
results and increase the electronic heating further. A method which goes beyond white noise and includes 
the change in electron-phonon coupling when the system is far from the equilibrium positions, e.g. close 
to bond-breaking, remains a challenge for the future. 
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V. CONCLUSION 



We have developed a semi-classical Langevin equation approach, which can be used to explore current- 
induced atomic dynamics and instabilities in molecular conductors. The Langevin approach can be 
solved in the harmonic approximation to obtain eigenmodes and their excitation in the presence of 
current, as well as used for molecular dynamics simulations based on the full anharmonic potential. 
Our simple, approximate MD simulation indicates that anharmonic couplings play an important role 
for the energy redistribution and effective heat dissipation to the electrode reservoirs. However, the 
MD is computationally very demanding beyond simplified model electronic structures and interatomic 
potentials and further developments are necessary. We have used carbon chain systems both to illustrate 
the Langevin approach, and in order to high-light how graphene might offer a unique test-bed for research 
into current-induced dynamic effects. Especially, it is straight forward to employ a gate potential to 
the gate electrode, and thereby obtain independent control of current and voltage bias in the system. 
Furthermore, atomic scale resolution can be obtained in electron microscopes in the presence of current, 
and Raman spectroscopy can give insights into the excitation and effective temperature originating from 
the electronic current [5TJrl52j. Our results for the simplified carbon chain systems indicates that it may 
be possible to tune the current-induced instabilities in the atomic dynamics with gate and bias voltages 
in the experimentally relevant range. 
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